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1  Foreword 

This  report  regards  in  detail  the  research  carried  out  under  AFOSR  Grant  FA9550-06-1- 
0425,  “Discontinuous  Galerkin  for  Diffusion,”  in  the  final  project  period  of  9  months  (June 
1,  2007  -  February  29,  2008),  and  then  reviews  the  entire  yield  of  the  project,  thus  serving 
as  the  Final  Report.. 

The  original  award  period  for  the  project  ran  from  1  June  2006  through  30  November 
2007;  a  3-month  no-cost  extension  was  requested  by  the  PI  for  medical  reasons,  and  granted.. 
Thus,  in  the  second  year  of  the  project  the  grant  only  covered  the  9  months  from  1  June 
2007  through  29  February  2008. 

We  bring  into  mind  that,  as  reported  in  the  first  annual  report  covering  the  period  from  1 
June  2006  through  31  May  2007,  this  project  on  the  recovery-based  Discontinuous  Galerkin 
(RDG)  method  for  diffusion  operators  has  been  rich  in  early  successes.  Most  of  the  tasks 
listed  in  the  2005  proposal  were  accomplished  in  the  first  year;  in  addition,  an  unexpected 
fundamental  concept,  the  “recovery  basis,”  was  developed. 

2  Accomplishments/new  findings  in  the  second  year 

In  the  second  year  the  effort  addressed  the  following  issues: 

1.  Making  RDG  suitable  for  multidimensional  applications.  This  meant  coming  up  with 
the  proper  polynomial  bases  for  recovery  at  the  interface  between  arbitrary  triangular 
or  tetrahedral  cells.  The  2-D  rule  is  described  in  an  AIAA  paper  by  Van  Leer,  Lo  and 
Van  Raalte  [1];  the  3-D  rule  is  a  nontrivial  extension  and  still  awaits  publication. 

2.  Increasing  its  efficiency  by  constructing  the  recovery  basis  once  as  pre-processing  step 
at  the  start  of  an  evolutionary  or  steady-state  calculation.  At  the  basis  lies  the  Fun¬ 
damental  Theorem  of  Recovery  by  Van  Raalte  and  Van  Leer  [2],  which  says  that  the 
expansion  on  two  neighboring  elements  of  the  discontinuous  discrete  solution  in  terms 
of  the  original  discontinuous  basis  functions  is  identical  in  coefficients  to  the  expansion 
of  the  recovered  smooth  solution  in  terms  of  the  recovered  smooth  basis  functions  (the 
recovery  basis). 
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3.  The  combination  of  RDG  with  upwind  DG  for  advection  (URDG).  This  is  accomplished 
by  using  the  representation  of  the  solution  in  terms  of  the  original  discontinuous  basis 
when  computing  inviscid  fluxes  (with  a  Riemann  solver),  while  using  the  representation 
in  terms  of  the  recovery  basis  to  compute  the  diffusive  fluxes.  Numerical  tests  were 
carried  out  based  on  the  linear  advection-diffusion  equation  and  on  Burgers’  equation. 

4.  Writing  the  general  1-D  form  of  RDG  as  a  penalty  method.  After  rewriting  the  RDG 
formulas  for  a  number  of  polynomial  orders  as  a  classical  method  with  additional 
penalty  terms  with  unique  coefficient  values,  the  general  form  of  the  expansion  became 
evident.  The  penalty  terms  are  different  for  odd  and  even  polynomial-space  degrees  k ; 
they  contain  the  product  of  the  jump  [i/,(A]  of  the  kth  derivative  of  the  solution  at  the 
interface  with  the  jump  of  either  the  test  function  <f>  (for  even  k)  or  its  derivative  (j)x 
(for  odd  k). 

5.  Dissemination  via  conferences  and  publications.  The  results  of  the  project  were  pre¬ 
sented  at  the  2007  AIAA  CFD  Conference  in  Miami  and  the  2007  International  Confer¬ 
ence  on  Spectral  and  Higher-Order  Methods  (ICOSAHOM  08)  in  Beijing.  The  Miami 
presentation  is  contained  in  an  AIAA  paper  [1],  the  Beijing  presentation  is  to  appear 
in  a  special  issue  of  Communications  in  Computational  Physics[ 2]  in  2008. 

6.  Technology  transfer.  By  leveraging  interactions  with  WPAFB  in  the  framework  of 
the  Michigan/Air  Force/Boeing  Collaborative  Center  for  Aero  Sciences,  we  were  able 
to  get  in  contact  with  a  representative  of  HyperComp  (Westlake  Village,  CA),  and 
formulate  an  STTR  announcement.  This  was  adopted  and  issued  by  the  Computa¬ 
tional  Mathematcs  program  of  AFOSR.  (NB:  The  ensuing  HyperComp/UMich-alliance 
STTR  proposal  was  indeed  succsessful.) 

3  Accomplishments  in  the  entire  grant  period 

Over  the  total  grant  period  the  RDG  method  evolved  from  a  promising  one-dimensional 
Discontinuous  Galerkin  discretization  technique  for  diffusion  terms  with  superior  accuracy 
and  good  stability  to  an  efficient  multidimensional  methodology  with  a  solid  theoretical 
underpinning,  and  ready  for  transfer  to  Air  Force/industrial  use.  A  comprehensive  list  list 
of  accomplishments  is  given  below. 

(A)  Accuracy  in  one  dimension.  Consider  a  uniform  spatial  mesh  of  open  elements  ftj  = 
]jh,  (j  +  1 ) A. [ .  Assume  the  approximate  solution  u  is  represented  in  each  element  on  a 
hierarchical  polynomial  basis  of  degree  k.  If  we  want  to  compute  the  diffusive  fluxes  at  the 
interface  between  VLe  and  f2e+1,  the  recovery  principle  says  we  must  recover  a  smooth  solution 
approximation  /  on  (J  Qe+1 ;  this  function  lies  in  a  polynomial  space  of  degree  2k  + 1.  The 
discontinuous  solution  on  f2e  and  De+i  is  the  L2  projection  of  /. 

The  degree  of  /  suggests  that  the  spatial  order  of  accuracy  of  the  resulting  DG  operator 
will  be  2k  +  2;  eigenvalue  analysis  and  numerical  experiments  show  that  the  order  is  higher. 
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Specifically,  we  found  that  the  accuracy  of  the  cell  average  of  the  solution  is  of  the  order 
3k  +  2  for  k  even,  and  3k  +  1  for  k  odd. 

When  solving  a  1-D  Poisson  problem  in  the  piecewise  polynomial  space  with  k  >  2,  the 
coefficients  of  the  polynomials  of  order  <  k  —  2  are  obtained  exactly,  i.  e.,  within  round¬ 
off  [1], 

(B)  Accuracy  in  two  dimensions.  The  property  of  (k  —  2)-exaetness  is  lost  in  two  and 
more  dimensions.  We  verified  numerically,  by  solving  a  2-D  Poisson  problem  on  a  uniform 
rectangular  grid,  that  the  order  of  accuracy  is  2k  +  2.  It  has  been  communicated  to  us  by 
Hung  Huynh  (NASA  Glenn  RC)  that  with  the  use  of  a  tensor-product  basis,  rather  than 
the  minimal  complete  basis  of  order  k,  the  2-D  method  inherits  the  order  of  accuracy  of  the 
1-D  method.  We  have  not  yet  confirmed  this. 

( C)  Boundary  treatment  and  accuracy.  At  the  boundary  of  a  1-D  domain  there  are  only 
k  +  2  data  available  for  recovery,  if  one  just  includes  the  cell  on  the  boundary.  The  resulting 
reduction  in  accuracy  of  the  boundary  fluxes  contaminates  the  accuracy  in  the  interior  of  the 
domain;  this  has  been  confirmed  in  1-D  numerical  experiments  The  loss  can  be  avoided  by 
allowing  the  missing  k  data  to  come  from  the  next  cell  inward,  starting  with  the  lowest-degree 
information  [1], 

This  approach  also  works  on  multidimensional  grids,  for  recovery  along  the  normal  to 
the  domain  boundary,  as  has  been  confirmed  in  2-D  numerical  experiments  on  a  rectangular 
grid  [1] .  At  the  end  of  the  grant  period  we  started  experimenting  with  the  extension  of  this 
technique  to  triangular  grids. 

(D)  Two-dimensional  recovery.  When  applying  the  recovery  principle  in  two  dimensions, 
it  is  beneficial  to  switch  to  a  coordinate  system  aligned  with  the  interface  across  which  the 
diffusive  fluxes  must  be  calculated.  In  the  ^-direction,  normal  to  the  interface,  the  recovery 
problem  resembles  1-D  recovery.  If  in  the  neighboring  elements  the  solution  is  represented  by 
a  complete  2-D  form  of  degree  k ,  the  ^-dependence  of  the  recovered  smooth  solution  will  be 
of  the  degree  2 k  +  1.  In  the  /"/-direction,  along  the  interface,  there  is  no  increase  in  accuracy: 
the  degree  remains  k  [1]. 

(E)  Stability.  A  general  stability  proof  for  recovery-based  spatial  DG  operators  is  still 
missing,  but  numerical  evaluations  of  eigenvalues  using  arbitrarily  high-precision  arithmetics 
have  confirmed  stability  up  to  A:  =  6.  In  two  dimensions  stability  has  been  shown  on 
rectangular  and  right-triangular  grids  for  k  =  1.  If,  on  a  rectangular  grid,  tensor-product 
basis  functions  are  chosen  (we  have  not  experimented  with  this),  1-D  stability  implies  2-D 
and  3-D  stability. 

(F)  The  recovery  basis.  For  convenience  consider  again  the  1-D  case.  When  applying  the 
recovery  principle  to  each  of  the  original  basis  functions  on  h2e  (J  De+1  (these  are  polynomial 
on  one  element  and  zero  on  the  other  one),  smooth  functions  result  that  may  be  used  as 
basis  functions  in  which  to  express  any  recovered  smooth  solution  on  Oe  (J  Oe+i .  To  make 
the  recovery  polynomial  basis  unique  we  require: 
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1.  the  discontinuous  basis  functions  fa  are  orthonormal; 

2.  the  recovered  basis  functions  ipj  are  orthonormal  with  respect  to  the  discontinuous 
basis  functions,  in  the  following  sense: 


[_J  ^e+1 


faifjdx 


1,  i  =  l,...,2k  +  2,  j  =  l,...,2k  +  2,  i  =  j] 
0,  i  =  l,...,2k  =  2,  1  =  1, ...,2K  =  2,  i^j. 


With  the  above  normalizations,  the  expansion  of  /  in  terms  of  the  recovery  basis  on  u  f^e+1 
is  identical  to  the  expansion  of  u  in  terms  of  the  discontinuous  basis  [3,  2]: 

2k+l  2k+\ 

u(x)  =  ^2  ai^i(x)  =>  f(x )  =  ^2ai  (^*0*0,  X  e  Oe(JOe+i.  (1) 

1=1  1=1 

While  the  discontinuous  basis  figuring  above  is  hierarchical,  with  maximum  degree  k,  all 
recovered  basis  functions  are  polynomials  of  degree  2k+l.  Examples  of  both  bases  are  found 
in  Figures  1  and  2. 

When  switching  from  one  discontinuous  basis  to  another,  the  transformation  between  the 
bases  applies  equally  to  the  recovery  bases. 

(G)  Flux  formulas  in  one  dimension.  For  the  computation  of  the  diffusive  fluxes  the  interface 


Figure  1:  The  discontinuous  and  recovery  Figure  2:  The  piecewise  constant  function 
basis  for  k  —  1  and  cells  h20  and  with  and  its  counterpart  in  the  recovery  basis 
h  —  1.  for  k  =  0  to  3. 

values  of  /  and  fx  are  needed.  Since  /  is  expressed  in  the  recovery  basis  functions,  and  these 
span  the  polynomial  space  of  degree  2k+l  on  [J  f2e+i,  /  is  independent  of  the  discontinuous 
basis  used  in  the  elements.  It  then  follows  that  the  diffusive  fluxes  are  also  independent  of 
the  basis  used;  only  the  degree  k  matters.  The  fluxes  can  be  uniquely  expressed  in  terms  of 
jumps  and  averages  of  the  discontinuous  solution  and  its  derivatives  at  the  cell  interface. 
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The  fundamental  theorem  of  recovery  makes  it  attractive  to  compute  and  store  the  re¬ 
covery  basis  for  each  interface  once  and  for  all  at  the  start  of  a  numerical  experiment,  if 
on  a  fixed  grid.  To  compute  a  diffusion  flux  at  an  interface  no  recovery  process  is  needed; 
one  uses  the  recovery  basis  in  combination  with  the  coefficients  of  the  discontinuous  solution 
representation  on  the  neighboring  elements.  This  makes  RDG  computationally  efficient.  At 
the  end  of  the  grant  period  we  were  in  the  process  of  implementing  RDG  in  the  above  way 
on  an  unstructured  triangular  grid. 

(H)  Bilinear  forms.  The  fluxes  go  directly  into  bilinear  forms.  For  k  =  1  the  form  looks 
like  the  symmetric  DG  method,  augmented  with  Arnold’s  interior  penalty  term  and  a  new 
penalty- like  term  based  on  [^x][«j;].  For  k  =  2  another  penalty  term  appears,  based  on 
[(f)]  [uxx\  i  Continuing  this  process  for  some  higher  values  of  k ,  the  general  form  of  the  expan¬ 
sion  became  evident.  The  penalty  terms  are  different  for  odd  and  even  polynomial-space 
degrees  k ;  they  contain  the  product  of  the  jump  [m^]  of  the  kth  derivative  of  the  solution  at 
the  interface  with  the  jump  of  either  the  test  function  <f>  (for  even  k)  or  its  derivative  4>x  (for 
odd  k). 

The  asymmetry  of  the  penalty  terms  starting  with  k  —  2  leads  to  the  appearance  of 
complex  eigenvalues  with  negative  real  part  in  the  spatial  DG  operator,  rather  than  purely 
negative-real  eigenvalues.  This  does  not  appear  to  affect  the  stability  of  the  RDG  method. 

4  Numerical  illustrations 

The  following  figures  illustrate  some  of  the  findings  listed  in  Section  3. 

Figure  3,  taken  from  [1],  shows  the  mesh  convergence  of  solution  properties  for  a  1- 
D  Poisson  problem,  with  k  —  2.  The  error  norms  of  the  solution  coefficients  =  A 2u 
(undivided  second  derivative)  and  a-2  =  A u  (undivided  average  gradient)  indicate  6th-  and 
7th-order  accuracy;  for  cq  =  u  (cell  average)  we  would  expect  8th-order  accuracy,  but  the 
errors  are  machine- zero  for  any  mesh  width. 

The  property  of  ( k  —  2)-exactness,  demonstrated  above  in  one  dimension,  does  not  carry 
over  to  multiple  dimensions.  Figure  4  shows  properties  of  the  solution  to  a  2-D  Poisson 
problem  on  [0, 1]  x  [0, 1]  with  Dirichlet  boundary  data.  The  exact  solution  is 

U(x,y)  =  -{cos(27tx)  +  cos(2-7n/)  —  1}.  (2) 

The  numerical  solution  was  taken  to  be  piecewise  linear  ( k  =  1  in  both  directions)  on  square 
meshes,  and  obtained  by  marching  the  associated  unsteady  problem  till  convergence.  The 
boundary  was  accurately  treated  according  to  Section  3,  item  (C).  The  left  plot  shows  that 
the  method  is  4th-order  accurate,  even  when  measured  in  the  maximum  norm.  This  order  is 
consistent  with  2k  +  2. 

The  solution  coefficients  and  errors  shown  on  the  right  are  for  the  16  x  16  mesh.  It 
is  worth  mentioning  that  the  numerical  solution  exhibits  all  4  symmetry  axes  of  the  exact 
solution. 
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Figure  3:  Mesh  convergence  for  1-D  Poisson  problem;  k  —  2.  Left:  convergence  with  mesh 
size  of  error  norms  of  solution  coefficients  A 2u  and  A u.  Right:  convergence  of  u. 


Next,  our  most  advanced  result  on  stability.  On  a  structured  grid  of  right  triangles,  which 
lends  itself  to  Fourier  analysis,  we  have  numerically  obtained  the  DG  operator’s  eigenvalues 
for  k  =  1;  all  are  non-positive  real.  Stencil  and  spectrum  are  shown  in  Figure  5. 

Finally,  results  for  nonlinear  advection-diffusion,  viz.,  Burgers’  equation, 


dtu  +  dx 


i^dxxU, 


(3) 


are  just  becoming  available.  We  are  using  the  test  cases  of  Wang  [4],  in  order  to  facili¬ 
tate  comparison  with  his  Spectral  Volume  method.  Figure  6(b)  shows  the  steepening  of  a 
compression  wave  in  time  for  /t  =  1  according  to  the  exact  solution 

/  x  2sinh(x) 

u{x,t)  =  —  ,  x  - j- — —]  4 

cosh(x)  +  exp  (3  —  t) 

we  computed  this  solution  with  k  =  1  on  a  sequence  of  grids  and  plotted  the  convergence 
of  the  Li  error  with  mesh  width  for  the  final  time  t  —  4  (left).  The  scheme  appears  to  be 
4th-order  accurate,  which  is  surprising  since  upwind  DG  for  advection  is  know  to  be  only  of 
order  2k  +  1,  that  is,  3.  Apparently  the  diffusion  error  dominates  for  the  meshes  chosen,  as 
the  diffusion  constraint  on  the  time  step  drives  the  CFL  number  to  smaller  values  on  finer 
meshes. 
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Figure  4:  Mesh  convergence  for  2-D  Poisson  problem;  p  —  1.  Left:  convergence  with  mesh 
size  of  error  norms  of  the  mesh  average  a\  =  u.  Right:  computed  solution  coefficients 
ai  (average  value),  a-2  (average  x-derivative)  and  <23  (average  y-derivative);  exact  solution 
coefficients;  coefficient  errors. 
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(a)  Exact  solution  at  different  times. 


(b)  Grid-convergence  study  at  t  =  4. 


Figure  6:  Steepening  of  a  compression  wave  according  to  Burgers’  equation,  used  as  a  test 
case  for  the  DG  method  with  k  —  1.  The  time  integrator  was  a  5-stage  4th-order  Runge-Kutta 
method. 
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